On the first Sonine correction for granular gases 
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We consider the velocity distribution for a granular gas of inelastic hard spheres described by 
the Boltzmann equation. We investigate both the free of forcing case and a system heated by a 
stochastic force. We propose a new method to compute the first correction to Gaussian behavior in a 
Sonine polynomial expansion quantified by the fourth cumulant 02- Our expressions are compared to 
previous results and to those obtained through the numerical solution of the Boltzmann equation. It 
is numerically shown that our method yields very accurate results for small velocities of the rescaled 
distribution. We finally discuss the ambiguities inherent to a linear approximation method in 02. 
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Most theories of rapid granular flows consider a granular gas as an assembly of inelastic hard spheres and assume 
uncorrelated binary collisions described by the Boltzmann equation, with a possible Enskog correction to account for 
excluded volume effects [llllllllllll0,llll[l3,lil|ll. The deviations from the Maxwellian velocity distribution 
may be accounted for by an expansion in Sonine polynomials, and it is often sufficient to retain only the leading term 
in this expansion, quantified by 02, the fourth cumulant of the velocity distribution 2^, 5, 8, 13, 14J. The purpose of 
this paper is twofold: first, we present a novel route to compute 02, directly inspired from a method that has been 
recently proposed to compute with accuracy the decay exponents and non Maxwellian features of gas sub ject ed to 
ballistic annihilation dynamics 0,^^ (where particles undergoing free fiight motion disappear upon contact [TtL IT^ ) . 
In essence, this method considers the limit of vanishing velocities of the Boltzmann equation, and deduces 02 from 
moments of the velocity distribution that are a priori of lower order than those involved in the standard derivation 
0, 13 . We may consequently expect a better precision from this alternative approach, that is analytically simpler 
to work out. We also know that the velocity distribution is non Gaussian at high energies so that extracting 

the relevant kinetic information from the behavior at vanishing velocities seems a promising route. The second goal 
of this article is to discuss the ambiguities -common to both approaches- encountered performing computations up 
to linear order in 02, neglecting not only higher order Sonine contributions but also terms in 02, k — 2,3. Such an 
ambiguity has first been mentioned by Montanero and Santos 8] . 

Within the framework of the Boltzmann equation, the one-particle velocity distribution function f{v;t) for a 
homogeneous system free of forcing obeys the relation 



dtf{^i;t)=I{f,f), 



where the collision integral reads 



I{f,.f) = <y 
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In Eq. (T is the diameter of the particles, 6 the Heaviside distribution, V12 = vi — V2 the relative velocity of two 
particles, V12 = Vi2/?;i2, W12 = |vi2|) and a a unit vector joining the centers of the grains. The space dimension is d. 
The precollisional velocities v** and the postcoUisional ones are related by 
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with a £ [0,1] the restitution coefficient. If energy is supplied to the system, an additional forcing term is present 
in Eq. Q 8], but the general arguments and method presented below remain valid. To be more specific, we shall 
also consider the situation where the system is driven into a non equilibrium steady state by a random force acting 
on the particles 0, 0, Q ■ With this energy feeding mechanism, coined "stochastic thermostat" , the Fokker-Planck 
term ■Co^v/ should be added to the right-hand side (r.h.s.) of Eq. (QJ where ^0 is related to the amplitude of the 
random force acting on the grains. 

We are searching for an isotropic scaling solution /(c) of Eq. The requirement of a time independent behavior 
with respect to the typical velocity vo{t) = ^J2l^v^)jd imposes that 
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where the rescaled velocity is given by c = v/vo{t) and the angular brackets (•) denote the average over /(v; t). The 
presence of the density n on the r.h.s. of Eq. |0J ensures that / dcf{c) = 1 and / dcc'^f(c)= d/2. This scaling 
function describing the homogeneous cooling state satisfies the time- independent equation 0, H, |^ 



^(rf + ci^)/(ci) =/(/,/), 



where 



and 



dcic^/(/,/), 



lifJ) 



dc2 / da 9{a- ■Ci2){^ ■ c 



12) 
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It is useful to consider the hierarchy of moment equations obtained by integrating Eq. Q over ci with weight 4] 

(8) 



The solution of Eq. Q is non-Gaussian in several respects. The high energy tail is overpopulated compared to the 
Maxwellian ^ , a generic although not systematic feature for granular gases (a particular heating mechanism leading 
to an under-population at large velocities has been studied in j^). Deviation from Gaussian behavior may also be 
observed at thermal scale or near the velocity origin. To study the latter correction, it is convenient to resort to a 
Sonine expansion for the distribution function /(c) [3 



/(c)=M(c) 



(9) 



where M.{c) — tt exp(— c^) is the Maxwelhan, and Si{c^) the Sonine polynomials (that may be found in W\; the 
first few are recalled in Q). Due to the constraint (c^) = d/2 the first correction ai vanishes 0|, and for our purposes 
it is sufficient to know 6*2(3;) = /2 — {d + 2)x/2 + d{d + 2)/8. From Eq. Q and making use of the orthogonality of 
the Sonine polynomials with respect to a Gaussian measure, one may relate the coefficient 02 to the kurtosis of the 
velocity distribution 



d{d + 2) 



(02 + 1), 



so that, upon taking p = 4 in Eq. ©, we get 



^4 = (rf + 2)(1 + 02)^*2 



(10) 
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In the following analysis, we will only retain the first correction in the expansion (Q: / = M{\ + a2>5'2). Computing 
/i2 and to linear order in 02 with this functional ansatz [and further linearizing Eq. I|ll|l ]. one deduces a2 0,0- 
This approach is nonperturbative in the restitution coefficient. However, since the high energy tail of A^(l -I- 025*2) 
is very distinct from that of the exact solution of Eq. (jsj , computing 02 from relation ((SJ with p > 4 is expected to 
give a poor estimate, all the worse as p increases. With this in mind, it appears that the limit of vanishing velocity 
of the rescaled Boltzmann equation contains an interesting piece of information: 



^i2f{0)= lim /(/,/). 



(12) 



The main steps to compute this limit are given in appendix. Up to a geometrical prefactor, the loss term of lim / on 
the r.h.s. reads /(0)(ci) and is thus of lower order than the quantities appearing in 111|) . Working at linear order 
in 02, one may therefore expect to achieve a better accuracy when computing the various terms (except may be the 
gain term) appearing in (|12|l than in (|ll|l . In the context of ballistic annihilation, a related remark lead to analytical 
predictions for the decay exponents of the dynamics and non-Gaussian features of the velocity statistics, in excellent 
agreement with the numerical simulations |l5l Il6j | . In the present situation, the gain term of / in (|12|) cannot be 
written as a collisional moment, so that the situation is less clear and deserves some investigation. We propose to 
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FIG. 1: Comparison of the correction 02(0;) for the free cooling in two dimensions obtained in Q, with Eq. H13^ . The crosses 
correspond to the "exact" result, obtained by solving the Boltzmann equation with the DSMC method, for 10® particles and 
approximately 500 collisions for each particle. The inset is a zoom in the region of the smallest root of the fourth cumulant. 

compare the value of 02 following this route to the standard one of Refs. Ulilil- Evaluating ^ at first order in 
02, we obtain: 



02 = 



4(a2 + l)^(a^-l) [72(a^ + l)-2] 
A{a,d) 



where 



(13) 



A{a, d)^5 + d{2-d) + 8a{a^ + l)(rf - 1) - ^^(23 - 6d + d^) + a*(3 + 6rf + d^) 

+ a^{-l + 2d + d^) - V2{a^ + ifia^ - 1)(3 + Ad + 2d^)/A. (14) 

In Fig.^ we compare this result with the analytical expression of van Noijc and Ernst 4J. We also display the fourth 
cumulant 02 obtained by Monte Carlo simulations from the numerical solution of the nonlinear Boltzmann equation 
(so called DSMC technique 113, |^). Our expression appears more accurate at small inelasticity, but less satisfying 
close to elastic behavior. The smallest root of 02 = obtained with Eq. H13|) is a 

root differs from the value a** — l/\/2 ~ 0.707. . . obtained upon solving (both a* and a** do not depend on 
space dimension d). The inset shows that the exact root is located in the interval ]a*, and seems closer to a** . 

In order to understand the discrepancy close to the elastic limit shown in Fig.^ it is useful to study the first Sonine 
correction f{ci)/A4{ci) = 1 + a2S2{cf). The result for a = 0.8 where our method seems to be the less accurate is 
shown in Fig.|21 and in Fig.Ofor a = 0.5. 

In spite of the imprecision of our analytical expression for 02 seen in Fig.Q] Fig. [21 shows that the limit method is 
very accurate for small velocities, but turns to quickly become more imprecise for bigger velocities. This suggests that 
computing the fourth cumulant from the limit of vanishing velocities gives more weight to this region which leads to a 
better behavior of the Sonine expansion for small velocities. On the other hand, the traditional route yields a global 
interpolation for all velocities. The good precision of our result for small velocities and the lower accuracy for higher 
velocities is confirmed in Fig. |31 Exploiting the above qualitative interpretation of the limit method, we expect to 
archieve a good accuracy using Eq. H13|) in order to find the first moment [T^ : 



(H> = f (1 
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Indeed, we suppose that the function 02 obtained from the limit method gives a precise description of the rescaled 
velocity distribution for small velocities. Thus our 02 is likely to describe more accurately a low order velocity moment 
than a high order one. This is confirmed by Fig. ^ 



4 




FIG. 2: Plot of f{c,)/M{ci) for a = 0.8. The curve labelled "Eq. (13)" and "Noije/Ernst" correspond to 1 + 02^2 where 02 
is given respectively by Eq. 1131 and by the Sonine correction obtained by Noije and Ernst following the traditional route A\. 
"DSMC" refers to the full distribution obtained from the solution of the Boltzmann equation (using 10^ particles and averaging 
over 300 independent samples). 




FIG. 3: Same as Fig.Elfor a = 0.5. 



As emphasized by Montanero and Santos [Sj , a certain degree of ambiguity is present when evaluating an identity 
such as (|llf) or H12|) to first order in a2. According to the way we rearrange the terms /i4, /Lt2, and {d + 2)(1 + 02) in 
say Eq. (|ll|l and subsequently apply a Taylor series expansion in 02, we obtain different predictions for 02(0). For 
instance van Noije and Ernst did expand the relation (|II|) 4], whereas Montanero and Santos also considered other 
possibilities such as /i4//i2 = (d + 2)(1 + 02) (this leads to a result which turns out to be fairly close to the one in 
U) and also /i4/(l + 02) = [d + 2)/i2. For small a in the latter case, the resulting 02 turns out to be 20% lower 
than the previous ones, and very close to the exact (within Boltzmann's equation framework) numerical results, for 
all the values of the restitution coefficient @. We push further this remark and show in Fig. the eight simplest 
different possible functions 02(0) obtained upon rearranging the terms of Eq. Hll|) and expanding the result to first 
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FIG. 4: First rescaled velocity moment (|c|) as a function of the restitution coefficient. DSMC is done for 10^ particles and 
approximately 500 collisions for each particle. 
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FIG. 5: The eight possible fourth cumulant 02 obtained from Eq. (IIH . corresponding to the two-dimensional homogeneous free 
cooling. We define rj = {d + 2){l + a2), then rewrite the equation jj,4 — ri^2 according to the eight possible different combinations 
mentioned in the legend, before doing the linear Taylor expansion around 02 = 0. The first curve is the plot of the function 
02 obtained by van Noije and Ernst Q], whereas the second one - obtained by Montanero and Santos igj- is very close to the 
exact results shown by crosses. 



order in 02. A similar ambiguity is present making use of Eq. 1)12(1 . The corresponding eight different possibilities are 
plotted in Fig. It appears that the envelope of the curves following from this method is less spread than within 
the "traditional" route, by a factor of approximately 2. We thus achieve a better accuracy at small a. 

The dispersion of the curves in Figs. |S1 and illustrates the nonvalidity of the linearization approximation at small 
a. However -and concentrating on Fig. (S]- it appears that all curves do not have the same status. Brilliantov and 
Poschel have indeed solved analytically the full nonlinear problem [i.e. working again with the distribution function 
/ = (1 + 0282) but keeping nonlinear terms in 02] , and obtained results that are very close to those of Noije/Ernst, 
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FIG. 6: Same as Fig. making use of Eq. 112|l instead of Ijlljl to compute the first Sonine correction. In the legend, I denotes 
lim/ and fo = 7(0). 



except for a < 0.2 where they found shghtly larger fourth cumulants Their result is therefore farther away 
from the exact one obtained by DSMC (see e.g. Fig. ^ where it appears than the Noije/Ernst expression already 
overestimates the exact curve). The difference between the DSMC results and those of Brilliantov/Poschel therefore 
illustrates the relevance of Sonine terms ai with j > 3 in expansion Q. However, some of the curves shown in Fig. 
lie close to the exact one, which means that it is possible to correct the deficiencies of truncating / at second Sonine 
order by an ad-hoc linearizing scheme. The agreement obtained is nevertheless incidental, and the corresponding 
analytical expression should be considered as a semi-empirical interpolation supported by numerical simulations. One 
should thus emphasize that the right way to compute 02 is to use its definition involving the fourth rescaled velocity 
cumulant of Eq. IjlOl) because this relation is not sensitive to higher order Sonine terms, nor to nonlinearities, even if 
this route doesn't give the most accurate description in the small velocity domain_(as seen from Figs. |21 and 

For completeness, we now briefly consider the stochastic thermostat situation 0,0,IS^3j where the counterpart 
of Eq. jSJl reads 



-gv^J(ci) = /(/,/). 

Considering again the limit ci — > and retaining only the first correction in the expansion (PJ, we get 



2^" 



02- 



(d + 2)(d + 4) 



- lini /(/,/). 



(16) 



(17) 



Given that the r.h.s. is already known from the free cooling calculation, it is straightforward to extend the previous 
results to the present case. As before, there are 8 possible ways to extract 02 from Eq. IjlTII working at linear 
order. The resulting expressions are displayed in Fig. |7| On the other hand, the moment method described in 
Refs. 0, makes use of the identity ^2{d + 2) = ^4, that is a direct consequence of Eq. (fTE|l . There are thus 4 
possible rearrangements leading to the different cumulants shown in the inset of Fig. For comparison, we have 
also implemented Monte Carlo simulations in the present heated situation (see the crosses in Fig. [T)). It is difficult 
to compare the dispersion of the curves with both methods (8 possibilities versus 4), since our approach makes use 
of Eq. H17|) which is of higher order in 02 than fi2{d + 2) — /i4, the starting point used in Refs. Our method 

appears here less accurate than for the free cooling, with again an underestimation of 02 at large a. 

In order to get free from the ambiguities inherent to a linear computation in 02, we have also solved the full nonlinear 
problem. The computation becomes cumbersome, and since Brilliantov and Poschel have already initiated this 
route in 3D for the homogeneous free cooling (thereby providing the calculation of 112 and /Z4), we will turn our 
attention to the 3D situation. First and for the sake of comparison, we have repeated the nonlinear derivation of Ref. 
for the stochastic thermostat. Second, we have computed the right-hand sides of Eqs. H12|) and H17|) without any 
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FIG. 7: The counterpart of Fig. |S|for the two dimensional stochastic thermostat. The inset shows the 4 possibilities associated 
with the method of Refs. 0,0]. The symbols show the results of DSMC simulations. 
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FIG. 8: Fourth cumulant in 3 dimensions for a force free system in the regime of homogeneous cooling. The curves correspond 
to the nonlinear solutions of Eqs. (IIH and 1121 (see text for details). The crosses correspond to the Monte Carlo results. The 
inset shows the same curves for the stochastic thermostat. 



linearization, from the form / = A^(l + a252). The left-hand sides only require the knowledge of /i2- For both free 
and forced situations, we subsequently obtain a polynomial equation of degree 3 for 02 from which we extract the 
physical root, the two others corresponding to unstable scaling solutions 0. The results are displayed in Fig. |H1 In 
particular, our approach again suffers from an underestimation of 02 for a > 0.5, already observed within the linear 
computation, and that is thus ascribable to Sonine terms of order 3 or higher. In this respect, it is surprising that 
these terms do not affect similarly the moment method of Ref. 9] in the same range of inelasticities. 

To sum up, using a new approach we obtain the first non-Gaussian correction 02 to the scaled velocity distribution. 
In view of the above results, we conclude that our approach constitutes an improvement over the previous procedures 
in the small velocity regime, and our analysis turns to be technically simpler to perform. We have also discussed the 
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ambiguities that arise 1) when restricting ourselves to second Sonine order, and 2) when a further hnearization of the 
various relevant relations is performed. It appears that an ad-hoc linearization scheme (point 2) may circumvent the 
limitations inherent to point 1. In any case, the computation of a non-Gaussian correction suffers from uncontrolled 
approximations that systematically need to be confronted against numerical simulations. 

We acknowledge useful discussions with A. Santos and A. Barrat. This work was partially supported by the Swiss 
National Science Foundation and the French "Centre National de la Recherche Scientifique" . 



APPENDIX: CALCULATION OF THE LIMIT ci ^ OF THE COLLISION TERM 

We define the loss term Ii and gain term Ig by 

Ii = - lim / dc2 \ da6{a ■Ci2){a ■Ci2)I{ci)J{c2) 



Ig = lim — 



dc2 / d^9{B-Zi2){^-ci2)f{cr)f{cT), 



(A.la) 
(A.lb) 



so that limci^o /(/, f) = Ii + Ig. 

Taking the limit ci ^ of the loss term yields the exact result 



where 



-/3i/(0)(c2) 



/3i = / da9{a ■C2){^ ■C2) = -rrr-— 



(A.2) 



(A.3) 



with r the gamma function. Within the framework of the Sonine expansion (jSJ, neglecting the coefficients a.^, i > 3, 
and making use of (16| 



a2\ r[(d+l)/2] 



Eq. (|A.2|I becomes 



SdMjO) 
20F 



1 + 02 



T{d/2) 
d{d + 2) 



^2 

1 - 



(A.4) 



(A.5) 



where 5'^ = 27r'^/^/r(d/2) is the surface of the d-dimensional sphere. 

Defining (3 — {1 + a) /{2a) > 1, the precoUisional rescaled velocities c** and postcoUisional ones are related by 



cj* = Ci - /3(C12 • 

C2* = C2 + /3(ci2 • a-)a-. 



The gain term ljA.lb|) thus becomes 

1 



^3 - ^2 



dc2 / da- e{a ■ C2)(ct • C2) / [/3(c2 • ct)ct] / [02 - /3(c2 • a)a] , 



(A.6a) 
(A.6b) 



(A.7) 



where the function / is isotropic. Performing the integration over C2 before that over a, we choose the x Cartesian 
coordinate as corresponding to the &■ direction. The velocity C2 is thus written C2 — c^^x + cj^, with Cx = (c2 
and c_L = C2 — CxS. G R"*"^. Eq. ljA.7(l becomes 



Ig ^ ^ I da I dc2 0{cx) Cx f {/3cxa) f {c2 - f3cxS-) 



Sd 
a- 



— I dcx 



dc^f{Pcx)f[Jcl+cl{l-(3y 



(A.8) 
(A.9) 
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Eq. ljA.9|l is an exact relation within Boltzmann's framework. Making use of the the Sonine expansion ^ where we 
retain only the first correction 02, Eq. HA.9I) becomes 



Sd 



Ap^Hi-P^U j rfe^e-i [i + a,S,{(i\l)\ {l + a,S, [d + cl{l - Pf]] . (A.IO) 



With the definition of the second Sonine polynomial 6*2(2;) = /2 — 2)x/2 + d{d + 2)/^, one sees that Eq. ljA.10|) 
may be expressed as a sum of products of the integrals 



J±(n) 
Jx{n) 



dcj_e "^cl, 

-1 



that may be computed using the general relation (a > 0) 



a(rf+n)/2 T{d/2) 



Tedious but technically simple calculations thus lead to 



SdMjO) 
2^ 



1 + 



+ a2Di{a, d) + a2-D2(a, d) 



(A.lla) 
(A.llb) 

{A.12) 
(A.13) 



where the final expressions Di{a,d) and Z?2(q:, c?) are too cumbersome to be given here. Finally, the limit ci — > of 
Eq. is given by the sums of Eqs. l|X5|) and lf03|) . 
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